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Abstract 


This report describes a stochastic collocation method to adequately handle a physically intrinsic uncertainty- 
in the variables of a numerical simulation. For instance, while the standard Galerkin approach to Polynomial 
Chaos requires multi-dimensional summations over the stochastic basis functions, the stochastic collocation 
method enables to collapse those summations to a one-dimensional summation only. This report furnishes the 
essential algorithmic details of the new stochastic collocation method and provides as a numerical example 
the solution of the Riemann problem with the stochastic collocation method used for the discretization of 
the stochastic parameters. 


1 Background 

With the continuous development and improvement of both CFD simulations reliability and accuracy, further 
developments toward physically relevant results raise the need to account for more realistic operating condi- 
tions. Past simulations used deterministic parameters such as for the boundary conditions, initial conditions, 
geometry of the problem, physical properties, etc. Meanwhile, the actual conditions are not known precisely 
and one must account for a certain level in uncertainty in the simulations. In a more general framework, it 
is indeed of critical importance to quantify uncertainty ( e.g . [9, 1, 17]) and establish a confidence interval 
for the simulation-based predictions or design ([3]). Unlike reliability analysis which deals with rare but 
catastrophic events, uncertainty- analysis is concerned with the probabilistic aspects of the simulations due 
to the stochastic nature of physical parameters, initial and boundary conditions. Some sources of uncertainty 
identified in numerical simulations include: 

• inaccuracy / indeterminacy of initial conditions (experimental data, variability of operating environ- 
ment), 

• bias in physical models (turbulence model), 

• approximations in mathematical model (due to linearization or assumptions that neglect some physical 
effects such as temperature and/or pressure dependence of relevant physical parameters), 

• uncertainty in describing the physical reality (errors in geometry, roughness, boundary conditions), 

• discretization errors (numerical diffusion, round-off, algorithmic errors due to incomplete iteration 
process, etc.). 

See [1] or [4] for a more complete review. 

One crucial issue is to accurately propagate uncertainty from, say, the boundary conditions to the bulk 
field and finally to the output of the simulation. Several techniques have been proposed in the past. Among 
these, the currently most commonly used if the so-called perturbation technique that basically involves ex- 
panding the variables of the problem in terms of Taylor series around their mean value (e.g. [11]). Though 
effective, this technique is limited to Gaussian or weakly non-Gaussian processes due to the difficulty in 
incorporating terms of order higher than two. Still, this technique allows for a quick estimation of the low 
order statistics. 

When one wants to get access to the full statistics, one can make use of the Monte Carlo technique. It 
is considered as an “exact” method for accounting for uncertainty in the sense that it does not require any 
approximations nor assumptions. The main advantages are that its convergence rate does not depend on the 
number of independent random variables and that its application is straightforward. However, it becomes 
intractable for large problems as it requires to earn- thousands of simulations which results in prohibitive 
CPU time. Indeed, the resulting approximated variance of a <r 2 -variance random process is a\ JC — o jn 
where n is the number of independent samples. Despite some techniques which improve the convergence 
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rate (Latin Hypercube, importance sampling, etc.), this last point prevents its use for most of the practical 
problems currently studied by CFD. 


Another possibly viable approach consists in deriving the inverse of the stochastic operator (that rep- 
resents the governing equations) by expanding it in a Neumann series. Meanwhile, the derivation of the 
Neumann series may however becomes difficult, if not impossible, for complex problems. 

An alternative, and more effective, approach is the Polynomial Chaos. Originating in the work of Wiener 
(1938) and mainly applied in the field of structural mechanics, this technique is now making its way in the 
fluid mechanics community ([15, 16, 12, 14]). The basic idea is to project the variables of the problem onto 
a stochastic space spanned by a set of complete orthogonal polynomials ^ that are functions of a random 
variable £(#), where 0 is a random event. The terms of the polynomials are functions of £(#) and are thus 
functionals. Many random variables can be used: for example, £(#) can be a Gaussian variable associated 
with Hermite polynomials. This is the original form of the W iener work and is called homogeneous chaos. 
Other expansions are possible: Laguerre polynomials with a Gamma distributed variable, Jacobi polynomials 
with Beta distribution, etc. Obviously, the convergence rate, and thus the number of terms required for a 
given accuracy, depends both on the random process to be approximated and on the random variable used. 
This aspect was illustrated by [15]. 


Using this approach, each variable of the problem to simulate is expanded as, say for the velocity u : 

oo 

u(x,t, 9) = ^(£( 0 )) . ( 1 ) 

I— 0 

For practical simulations, the series has to be truncated to a finite number of terms, hereafter denoted 
N PC . This framework remains valid also for partially, or non-, correlated variables. In that case, the variables 
of the problem are functions of several independent random variables and £ is now’ a vector. The general 
expression for the Hermite polynomial H is 




d p 

<9£„ . . . 


( 2 ) 


As discussed in [2], there is a one-to-one correspondence between functions H p (£ iu . . . ,£ in ) and v l' J (£(6 1 )). 
The general expression for the Hermite Chaos is finally given by 


N PC 

u(x,t,e)= *i(£(0», 

i=0 

with the number of terms Npc determined from 


(3) 


Npc + 1 


iP'pc “h Ppc ) • 
ftjDC • Ppc ■ 


(4) 


where p pc is the order of the expansion and n pc the dimensionality. It follows that N PC grows very quickly 
with the dimension and the order of the decomposition. As an example, for a second order, 2-D Hermite 
polynomial expansion, we get 


u(x,t,6) = «o(* ! O+wi(x,f)£i(0)+«2(a:,f)6(0)+W3(*,O (£j 1(0) - l)+«4(*,0? 1 (®)6W+« 5 (*,t) (£|(fl) ~ 1) 

, x (5) 

where £1 (0) and f 2 (0) are two independent random variables. 


In the case of the Hermite polynomials, the zero th order represents the mean value and the first order 
the Gaussian part while higher orders account for non-Gaussian contributions. The Hermite polynomials 
span a complete orthogonal set of basis functions in the stochastic space. The inner product is expressed as 
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( 6 ) 


</i(€) /a(€)> = / /i(€) /2(C) «>(€) rfC , 

./ —03 

where the weight function tw(£) is the n pc -D normal distribution: 


«■’(£) = 


y(2 tt)"*' 


e 2 * * 


(7) 


In the L 2 norm space, we thus get 


*j) = (*? ) *y - 


(8) 


where 5 is the Kronecker operator. 


2 Full Spectral Approach for Quasi 1-D Nozzle Flow 

To illustrate the Polynomial Chaos approach and highlight the computational burden of the multidimensional 
stochastic summations, the Polynomial Chaos technique is derived here for a quasi 1-D nozzle flow. Mathelin 
and Hussaini ([7]) have presented results for inviscid quasi 1-D nozzle flow solutions using Polynomial Chaos. 


The system of compressible Euler equations in conservative form is 

Q t + F x = S , 


(9) 


where 


( pA \ ( pu A 

puA F = I pu 2 A + PA 
\ pE A ) \ puE A + PuA 


0 

S=| P dA/dx 

0 


( 10 ) 


with p the density, u the velocity, A the cross-section area, P the static pressure and E the total energy. 


The pressure can be removed from the above equations making use of the following equation: 

P 1 2 
E ~ (7-1 )P + 2 U ' 


( 11 ) 


Equation (9) thus becomes 



The Euler equations are discretized in space using the spectral element method. This spatial discretization 
is consistent with the spectral expansion involved in the Polynomial Chaos technique. Equation (9) is 
projected onto a spectral and stochastic basis and the Galerkin technique is applied. The spectral space is 
spanned by Chebyshev polynomials on points (in local coordinates x <E [—1,1]). 


x n = cos (^P) Vn e [0; A r ] C N , (13) 

where N is the number of points within each element. The interpolants of variables, say u, in spectral space 
read as 

N 

u(x,t) = Y,u n (t)h n (x). (14) 

n=0 
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where 


h n {oo) 'y " _ _ T m (x n ) T m (x) 

iV Cf) Cry} 

m=0 n m 


(15) 


with T m the Chebyshev polynomials and c* such as 


Ci — 2 i€{0;N} 


(16) 


Ci = 1 i <E]0; N [ C N , 


(17) 


Similarly, in the stochastic space, 


Npo 


i(0,x,t) = ^ Ui(x,t) ^*(£(0)) . 


(18) 


After rearrangement, it finally follows: 


{ N pc Npc N 

in; Pi, n A j VjVt) (h n , h p ) 

t— 0 j=0 n=0 

Npc Npc Npc N N 

«=i inn Pi.n Uj,m Ak 'ifj ty) (h n; h m ,h p ) 

i=0 j = 0 0 n=0 m=0 

Npc *Vpc Npc N N 

,iiiii Pi,n A k 'h j ^ k ^ l) (h n ,/l m ,/lp) 

V ?=0 j= 0 fc=0 n=0 m= 0 


V {Z € [0; A t pc], P € [0; iV]} C N 


(19) 


/ 0 


s = 


Npc N pc N pc -V /V 


(i-^EEEE ^,m ^ <*<*;¥* *i> (h n ,h m ,h p ) 


dx 

z=0 j— 0 fc— 0 n=0 m=0 
iVpc A T pc -^Vpc Npc A r A T A f 

Y~ El El EE El El El El ft-" «J.m u *.o -Q- 1 <*» ** *r *1) (h n ,h m ,h 0 ,h p ) 

t=0 j= 0 k = 0 r=0 n=0 m— 0 o— 0 ^ 


( 20 ) 


4 


/ jV p f A' p c iV p c JV A r 

E E S L L p,n Uj ' m ( 9i ^ 

0 j=0 k-0 n— 0 m=0 

dh r 


OF 

dx 


q 4 . 

(hn,h m *h p ) 4- A k 




( u dh m , 

ar-^ 


NpC A r pc A " PC A r A r A’ 


y y y y y y y] pt,u ^jt, 0 (^, 


- , h Q , fap^ 4“ ^hri*. hm « • ^p^ ^ 


-0 j=0 k ~ 0 r=0 n= 0 m=0 o=0 

• 4 * (("a7 ,/lm,,lu ’ /lp ) + { hn '~dx 

N PC Npc N pc A r A 7 

(7 — 1) ^ ^ ^ ^ ^ Pt.n Ej,m i Vk ^/) 

i=0 j — 0 A-=0 n=0 m=0 

A, ((^pA-’O + fa 


(h n ,h m ,h u ,hp) ^ + 


+ 


I m=U 

dhm . \\ 

ln ’ dx ,hp )) 


/ j i , > d-4fc , 

{hmhm'hp) < 


A r p<7 Ape Ape A T pc A A r A r 

El £ £ E E E E " p, - n Uj - m Ek '° * k * r 

i=0 j=0 k = 0 r=0 n=0 m=0 o= 0 


(fan . ^I nv , /lo ? ^p) 4* 


^ ((i£ ,hm}h ° ,hp ) + { hn 'i^r' h " ,hp ) + ( hn - hm ^- hp )) 


*v - 1 


A r pc A r pc A’ pc A r pc A'pc A r N A" A 


E E E E E EEEE Uj - m Ufc -° Ur « ^ ** ^ *« *«> 


r=0 jf =0 k= 0 r — 0 5=0 n=0 m=0 o=0 q = 0 

( hin , j ho ? j hp ) ^ 4~ -4. A 


dx 


V 


( , , dho 

I fin , flm? 5 


/iq , hp J 4* ( ? hm s ho 


— s hm, ho< hqi hp'j 4- ^ > ho , j ^p^ 


+ 


ah, 

ax 


where 7 is the specific heat ratio (7 — 1.4 for diatomic gas). 

The scalar products in the spectral space are defined by, say for (h n , h m , h 0 ) 


( 21 ) 


(h n , h m) h Q ) — J h n (x) h m (x) h 0 (x) dx 

= J, V V V Ta ^ rb ^ m) Tc ^ [ l Ta(x)T b (x)T c (x)dx ■ 

AT 3 “ “! “ C a 4 C c C n C m Co 7- 1 


a=0 6=0 c=0 


( 22 ) 


The nonlinearity in the inviscid momentum flux results in a 9-dimensional summation in Eq . (21 ). Five of 
these sums are due to the Polynomial Chaos expansion, and would remain even for a simple finite difference 
spatial discretization. For the extension of this problem to turbulent, viscous flow (see Mathelin, et al [8]), 
the corresponding equation for the turbulent kinetic energy contains a 12-dimensional sum, with seven of 
these due to the Polynomial Chaos terms. 


For this problem all the nonlinearities are of polynomial type. For the Riemann problem, a fundamental 
component of many modern algorithms for solving Euler and Navier-Stokes equations, the nonlinearities are 
more complex. The appropriate approach to this for Polynomial Chaos methods is not obvious. 
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3 Stochastic Collocation Method 


The stochastic collocation method was developed to enable application to Spectral Discontinuous Galerkin 
Methods (for which solution of the Riemann problem is a necessity) as well as to reduce the cost of Polynomial 
Chaos methods for more classical discretizations, such as described above for quasi 1-D nozzle flow. 

Let a and b be two independent random variables. In the Polynomial Chaos method, the infinite series 
for each variable would be expressed as, say, for a, 


oc 

G = ^ ] a i^i (0 ) 
2 = 0 


(23) 


and the product ab would be 

30 OC 

<“=EE OibjVitfWjiO • (24) 

i=0 j=0 

(Here we drop any dependence upon x and t and concentrate on just the random variable £ without reference 
to the random event parameter 9.) 

The finite series for a is then 

jVpc 

« = Y. **<(0 • (25) 

2=0 

The usual Galerkin truncation yields for the expansion coefficients of the product c = ab 

Npo A t pc 

(ab) k = E E a t bj (V, 'I'j S'*) Vfc e [0; Af pc ] C N . (26) 

i=0 j- 0 

Therefore, computing the Npc T 1 coefficients of the product ab takes of order A r J» c operations. The 
complexity is even greater for cubic products. Moreover, for non-polynomial nonlinearities, such as those 
that arise in the Riemann problem, the Galerkin truncation is not obvious. 

For traditional deterministic problems, these difficulties are often treated by resorting to a collocation 
method instead of the Galerkin method. The usual spatial collocation approach consists of projecting the 
equations into the physical space (x). But a “physical space” corresponding to the random variable f has 
limited physical meaning, and the spatial scheme cannot be directly transformed. To deal with this problem, 
we define a new stochastic space whose properties are known. 

In the stochastic collocation method (SCM) one lets the Probability Density Function (PDF) of the 
random variable £ serve as the basis for the transformation between the physical random variable £ and 
its artificial stochastic space. We use a to denote a point in this artificial stochastic space. The range of 
the corresponding Cumulative Distribution Function (CDF) is [0,1], and provided that the underlying PDF 
is non-zero in the interior of its domain, the CDF is strictly monotonic and therefore this transformation 
is bijective. Denote the CDF of £ by -^(O* We prefer to map the CDF to the standard domain [-1,1] of 
orthogonal polynomials, and so we define 


m ) = 2^(0 - 1 . (27) 

Finally, denote the inverse of by G$. Let a be a point in the domain [-1,1] of G. This is the new, known 
stochastic space that is the foundation of the SCM. We have then 

a : = F*(0 = 2^(0 - 1 (28) 

and 

? = <?€(<*) . (29) 

The SCM utilizes collocation points — 0, . . . , N qy in [-1,1] . These have associated quadrature weights 
u’i, chosen to give the best approximation to inner products on [-1,1]. In the present work we make the 
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particular choice of the Gauss-Legendre points and weights. The collocation points in the stochastic space 
are associated with the points = 0,...,Aq, in random variable space according to (29). 

Let us return now to the example of evaluating the product of two random variables a and 6. First, we 
construct the appropriate mapping functions, F a and F i, (and also their inverses G a and Gt>), from the PDFs 
of a and b. Second, we evaluate a and b at the points a t , obtaining 

0-i = G a (ck i) i = 0, . . . , A q (30) 

hi = G b (ai) ?=0 N q . 

These functions are then approximated in a-space by their interpolating polynomials, as for example, 


% 

a A<? (a) = , 

j-o 


(31) 


where the functions h % are the classical Lagrange interpolants based on the N q T 1 collocation points, i.e., 
hj is a polynomial of degree N q and hj(cki) = &ij. 

The representation of the product c = ab associated with a Galerkin method would be 


N q N q 

< a Nq b Nq h k >= < hihjhk > . 

i=0 j=() 


(32) 


But in the SCM we resort to the quadrature rule to approximate this a s 

N q N q N q 

< a Nq b Nq h k >^ hi(ai)hj(ai)hk(ai) m 

2=0 i =0 j =0 

Nq N q A/q (33) 

) SuSjiSu 10 1 

1=0 i = 0 j=0 

= akbkWk ■ 


So, the representation of the product is just 

N q 

G c (a) = c Nq {a) = Y'(a k b k w k ) h k (a) . (34) 

fc=0 

The last step is to map this result back into the random variable space to obtain the approximate CDF of 
c = ab. Recall that (34) gives the representation in a-space of the random variable c. We now need to back 
out the CDF of c from the transformation rules (28) and (29). (In these equations we make the appropriate 
identification £ = c.) We have, finally, for the CDF , J- c , of c 

< 35) 

«i(c;V + >). 

recalling that F c is the inverse function of G c . The PDF of c is then readily obtained by differentiation of 
(35). 

In our present implementation of the SCM we use standard interpolation methods for the mappings (28) 
and (35). For the forward mapping we start with a uniform distribution of the £* (a, in our example) for 
representing the PDF of f . We base the interpolant for the inverse mapping on the collocation points a,, in 
which case we have 

G c (<^i) = QibiWi 2 = 0,..., A q . (36) 
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since h k (cti) = S ik . 

The complexity of this implementation of the SCM in one dimension is N q , regardless of the form of the 
nonlinearity. In contrast the Galerkin method already scales as Np C for a simple quadratic nonlinearity. 
The quadrature in the stochastic collocation method does introduce additional errors compared with those 
of the Galerkin method, but these are reduced as N q increases. This is usually not a major issue since the 
collocation scheme allows for dramatic speed-up compared to full Galerkin approach and a relatively high 
value for N q is affordable. 

As an example of this process, consider the cubic equation 

fW) = u 3 , (37) 

where u is a Gaussian random variable, with mean 1.0 and standard deviation of 1/(2 0T. Figure 1 shows 
the PDF of u, the a-representation of a, the a-representation of f(u) = u 3 , and the corresponding PDF of 
/(u). 


j \ 

i i 


^ ■ V 



(a) PDF of input function u. 


(b) Projection of the input function into the Q-space. 
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(c) Output function f(u) in the a-space 


(d) PDF of output function f(u). 


Figure 1: Illustration of stochastic collocation for f(u) = u 3 . 
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4 Application to Stochastic Riemann Solver 

To illustrate the applicability of the above method, the stochastic Riemann problem is considered here. The 
field variables are assumed to be discontinuous at a time t and this leads to the generation of shock wave, 
rarefaction and expansion fan. See Landau fc Lifshitz ([5]) or Liepmann & Roshko ([6]) for a review of 
the physical phenomena involved. To compute the resulting state for later times, one can make use of the 
algorithm proposed by Sod ([10]). It produces strongly non-linear equations including rational powers and 
conditional expressions. While straightforward in a deterministic space, the application of such an algorithm 
becomes very tricky with a stochastic approach. For example, two expressions are to be applied depending 
on the pressure ratio: if P*/P > 1 , then use equation (a) and if P* /P < 1, then use equation (b). But 
when both pressures are non-deterministic, what does P* /P > 1 or P* fP < 1 mean ? 

Using the collocation scheme, the pressure ratio now has a deterministic expression in the a-space and 
it is possible to determine the resulting state. Figures 2 through 4 show the propagation of the stochastic 
description of the state across the discontinuity, z.e., the stochastic collocation solution of the stochastic 
Riemann problem. In these computations N q = 200. 




(a) Left state. 


(b) Interface state. 


(c) Right state. 


Figure 2: Density PDF evolution throughout a discontinuity. 





(a) Left state. 


(b) Interface state. 


(c) Right state. 


Figure 3: Velocity PDF evolution throughout a discontinuity. 


Preliminary results from Mathelin, et al ([8]) indicate that substantial CPU savings have been achieved 
by the stochastic collocation method in Polynomial Chaos solutions of quasi 1-D nozzle flow. In their simplest 
case, Euler flow, the computational time was reduced by a factor of 4. The computational savings is strongly 
dependent upon the discretization parameters. More details are available in [8] . 
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(a) Left state. 


(b) Interface state. 


(c) Right state. 


Figure 4: Pressure PDF evolution throughout a discontinuity. 
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